Simulating Open Quantum Systems using the Simple Lie Group SU{4) 



Minghui Xu,' D. A. Tieri,' and M. J. Holland' 

^JILA, National Institute of Standards and Technology and Department of Physics, 
University of Colorado, Boulder, Colorado 80309-0440, USA 
(Dated: February 27, 2013) 

We show that open quantum systems of two-level atoms symmetrically coupled to a single-mode photon field 
can be efficiently simulated by applying a 5t/(4) group theory to quantum master equations. This is important 
since many foundational examples in quantum optics fall into this class. We demonstrate the method by finding 
exact solutions for many-atom open quantum systems such as lasing and steady state superradiance. 
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Although quantum mechanics actually describes the behav- 
ior of isolated and closed quantum systems, in practice most 
physical situations to which it is applied are open. The open 
nature is necessary to treat basic irreversible processes such 
as the energy transfer with a heat bath, the particle or charge 
exchange with a reservoir, and quantum measurements. Open 
quantum systems can be treated under the Born and Markov 
approximations by the quantum master equation in the Lind- 
blad form lU, which has been applied across many fields of 
physics, including quantum optics and quantum information 
science Hlst], atomic and molecular physics ||4i], solid state 
physics ystl, and optomechanics j^. 

Here we focus on a novel group-theoretic approach to find 
an efficient solution of the quantum master equation for an 
important class of quantum optical systems. Even though we 
focus on these systems, the methods we present could be more 
generally applied to other fields. We consider specifically the 
symmetric coupling of a single-mode cavity field to an ensem- 
ble of two-level atoms (analogous to pseudo-spin- 1/2 sys- 
tems or qubits). The Hamiltonian that describes this situation 
in the interaction picture is given by 



(1) 



where the first term is the free energy, with A being the detun- 
ing of the light field from the atomic transition, and the second 
term is the reversible atom-field coupling with strength Q. The 
photon annihilation operator is a, and crj^ and cr| - (crj)^ are 
Pauli operators for the jth spin-component. In the presence of 
decoherence, the full quantum evolution is described by the 
quantum master equation for the reduced density operator p: 

p^£p= ^[H,p]+KD[a]p 
in 

+ \7!D[CTJ] + wD[a)] + —D[cr'j]jp , (2) 

where D[d]p = (lOpO'' -d^6p-p& d)/2 denotes the Lind- 
blad superoperator We have introduced the decay rate k for 
the cavity, and population relaxation rates for the spin compo- 
nents y, w (for decay and pumping respectively) and dephas- 
ing rate l/(2r2). 



In general, for all but smallest values of A^, exact ana- 
lytic solutions to Eq. (|2| are intractable, and it is neces- 
sary to resort to numerical simulation approaches such as the 
quantum trajectory method fl\. Even so, treating more than 
about ten spins is difficult due to the exponentially increas- 
ing dimensionality of the underlying Hilbert space (scaling 
as 2^^). It should be emphasized that the reduced basis of fully- 
symmetrical Dicke states do not suffice to solve the quantum 
master equation even though the Lindblad operator treats the 
spin-components symmetrically. As will be shown below, all 
state vectors with mixed symmetry are in principle required . 

Recently, it was pointed out that this is not true in Liouville 
space since the Lindblad operators are invariant under SU{4) 
transformations JH]. This observation allows one to express all 
of the Lindblad operators in terms of generators of the SU(4) 
group. For this purpose, 18 superoperators 0+, O and 
where O e {Q, 2, Ai, N, H, 'V] are defined and can be shown 
to be linear combinations of the 15 Gell-Mann matrices of the 
5'C/(4) Lie algebra (see supplementary material). 

As a consequence, it is possible to construct a reduced basis 
for the density operator using a multiplet of the 5'C/(4) group. 
Transcribing notation from the four-flavor quark model — a 
model that has the same symmetry structure — the fundamen- 
tal representation is given by u - \1}{1\, d = |0){0|, s - |1)(0|, 
and c - |0)(1| (up, down, strange, and charm). Since the sym- 
metry type of the basis is preserved under the action of the 
SU(4-) generators, this leads to a tremendous reduction of the 
number of required basis states needed to provide an exact 
solution of the master equation. 

For the fully symmetric case, the basis is: 
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where iS denotes the symmetrizer and a + (3 + y + 5 - N. 
Note that only basis states with y - S - Q have non- 
vanishing trace. The three quantum numbers q,q^ and 0-3 
have ranges q - 0,1/2,..., N/2, q^ - -q,-q + l,...,q and 
<T3 = q-N/2,q-N/2 + I, ...,N/2 - ^, resulting in the dimen- 
sionality of the basis (N + l)iN + 2){N + 3)/6, i.e. of order 
A^^*. This tremendous reduction should be compared with the 
case in which one does not take advantage of the accessible 
subspace and instead uses the full dimensionality of the Li- 
ouville space given by 4^^. Note that all the superoperators 
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commute with three associated Casimir operators (analagous 
to the S ^ total spin operator for SU{2)), and this makes it pos- 
sible to uniquely label the symmetry of the invariant subspace 
of Liouville space in terms of the Casimir eigenvalues. 

In this paper, we apply the SU{4-) group theory to find ex- 
act solutions to the quantum master equation in general form. 
We show how to calculate the various basic observables of in- 
terest. We demonstrate that the density matrix in the SU{4-) 
basis representation can be precisely mapped to the collec- 
tive spin-angular-momentum representation \S,M} in Hilbert 
space, which enables us to efficiently diagonalize the density 
matrix. This allows us to provide complete information about 
the system, including functional properties of the density op- 
erator such as the purity and von Neumann entropy. 

In order to solve Eq. Q, we expand the density matrix as 

p= Z c-;;;_^p,,,„<,,|m)(«|, (4) 

where C™^', are complex coefficients, and |«) is the photon 
Fock state. The Lindblad operators can be written compactly: 

N 
N 

D[(rf] = 4M3 - 2<33 - 2E3 - . (5) 

The completeness of O+ - 3 and a implies that an arbitrary 
Hamiltonian can be expressed by them, e.g. from Eq. dTJ, 

—[H,p] = -2iAI.2p - in \a{M+ + N+)p + a\M- + NS}p\ 
in L J 

+ iO. [('ZY+ + 'y+)pfl'' + CZY- + ^-)pa\ . (6) 

Combining Eqs. ^ and ^ with the action rules of the SU{4-) 
and photon operators on the basis states (see supplementary 
material) gives a closed solution of Eq. In general, this can 
be solved analytically or numerically with standard methods. 

Having established the procedure for determining the time 
evolution of p, it is now important to describe how to calculate 
physical observables. We begin with the trace of the density 
operator, which is given by 

Tr[p] = J] C-'^3 „ = 1 , (7) 

which is an invariant during evolution to represent probability 
conservation. Average values (a) and (a ' a) are found analo- 
gously. For the spin-operators, we provide the following ex- 
amples up to quadratic order: 

(trf > = 2TT[Q3p]/N, 

<crf frf ') = (4Trml - - A^)/[A^(A^ - D], 

{crp^Tr[{M^ + N^)p]/N, 

((r)cr, ) = Tr[n^_(At- + N-)p - Q-p\im^ - 1)], 



where j + k. We do not provide j-k quantities since they fol- 
low easily from the basic algebra of the Pauli spin-operators. 

For coherence properties it is necessary to calculate prod- 
ucts of operators evaluated at different times. Of particular in- 
terest are the first-order and second-order correlations, which 
can be found by applying the quantum regression theorem: 

<Oi(f + T)d2{t)} = Tr [()ie^^[O2p(0]l , 

{ddt)di{t+T)d2(t+T)62(t)} = Tv[62e^'[d2Pit)di]di] , (9) 

where e-^^[p] is the time propagation from Eq. Q starting with 
the initial density matrix p. For example, in order to obtain the 
first-order correlation of Oi and O2, one takes 02p{t) as an 
initial condition, time evolves it for t according to Eq. (|2]i, ap- 
plies 0\, and computes the trace. A similar procedure follows 
for the second-order correlation. In this way, field quantities, 
<fl"'"(f + T)a(t)) and {a\t)a''(t + T)a(t + T)a{t)) are direcdy cal- 
culated. For spin-coherence, the required expressions are: 

N 

^ <trt(r + T)(Tl{t)) = Tr + N^)e^'[{M- + N^p{t)]\ , 

j.k=i 

N 

^ {cr^it)cr} it + T)cr,{t + T)cr,, (f)) = 

Tr ['V-iM- + N-)e-^^['y-{M- + N-)p]] . (10) 

Although at this point we have provided a theoretical 
framework that is complete and provides exact and efficient 
solutions to the general quantum master equation, it is of- 
ten inconvenient to work in the Pq,q,,cri representation of the 
density operator For example, it can be a nontrivial proce- 
dure to characterize the many -body spin-state in this repre- 
sentation by quantifying the degree of entanglement, which 
is derived from a functional (i.e. Tr[plog(p)]). For this rea- 
son, we illustrate now the procedure for efficiently projecting 
the density operator from the SU{4-) basis representation onto 
the usual representation of density matrices formed from the 
Hilbert space basis vectors. These Hilbert space basis vec- 
tors are specified by the angular momentum eigenket 15, M), 
where S = N/2, N/2 - 1, ...,(l/2or0) is the total spin and 
M = -S , -S + 1, . . . ,S is the spin-projection. Note that S also 
labels the symmetry of the states, e.g. S = N/2 corresponds 
to the fully symmetrical Dicke states. 

In order to illustrate how this projection is done, it is in- 
structive for us to first examine explicitly the = 2 case where 
the Hilbert space is 4 dimensional. Two spins form a symmet- 
ric triplet state and an antisymmetric singlet state, correspond- 
ing to total spin 5 = 1 and 5=0 respectively. In this case, 
the complete density matrix from Eq. Q for given m, n is 
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Notice that the resulting matrix is block diagonal in the S - I 
and 5=0 subspaces (a 3 x 3 block and a 1 x 1 block). In 
addition, the complex coefficients contributing to the matrix 
elementfor|5,M)(5,M'| all satisfy ^3+0-3 =Mandg'3-cr3 - 
M'. Finally, the trace is simply 2;^,=_i C^q-ic = 1- 

These results can be systematically extended to higher A^. 
For any A^, the density matrix is block diagonal in 5, with 
each block given by 

p7-Yu^'s:mm\S^M){S,M'1 (11) 

MM' 

where Ds^.m' ^re density matrix elements for the symmetry 
type S . There are ns ways for spins to construct the ba- 
sis for each 5, so that 2^(25 + l)«s - 2'^, i.e. the Hilbert 
space dimension. To find its , we note that \S, M) forms a basis 
of the (25 + l)-dimensional irreducible representation of the 
SU{2) group. Determining «s is accomplished with the help 
of the Young tableau of the SU(2) group, where one can obtain 
the number of equivalent representations iteratively. Fig. [Ha) 
shows the Young tableau for the N - A case. A correspond- 
ing tabular method for evaluating iis for any is shown in 
Fig. [nb)^ which contains about one half of the well-known 
Pascal's triangle. 



(a) 



S=2 (n,= l) 



[— [— I S=0(n,=2) 




FIG. 1. (a) Young tableau for determining the irreducible represen- 
tations contained in the direct product representation of the 5(7(2) 
group for A' = 4. The dimension for 5 = 2, 1,0 is 5,3, 1 re- 
spectively, and n, = 1,3,2, so the total Hilbert space dimension is 
5x1-1-3x3+1x2 = 2'* as expected, (b) "Pascal's triangle" to 
evaluate ns in an iterative way for any A' (the case considered in (a) 
is the fourth row down from the top), (c) Pyramid representation of 
the density operator in the |S,M) representation for A' = 4. Each 
layer of the pyramid represents the block matrix for each 5 . 



Being able to express the density operator in the \S,M) rep- 
resentation makes easy the calculation of functionals, such as 
the purity Tr[p^], or the von Neumann entropy 



S - -Tr(p Inp) = - ^ A.j In Aj, 



(12) 



where Aj are eigenvalues of p. The point is that, because the 
density matrix is block diagonal in the \S,M) representation, 
we do not need to diagonalize the whole density matrix, which 
would be a daunting task. Instead, we only need to diagonal- 
ize a series of LA^/2J + 1 blocks of dimension 2S + \. 

In the following, we demonstrate the method by solving 
many-atom open quantum systems such as lasing and steady 
state superradiance. We show the capability for finding exact 
solutions of large systems and are able to obtain full informa- 
tion about both the transient and steady-state density matrix. 

First, let us consider a single-mode laser consisting of an 
ensemble of two-level atoms coupled to an optical cavity, 
which can be modeled by the general quantum master equa- 
tion Eq. (|2]i fl. In this model O is the atom-cavity coupling, 
K is the laser output coupling loss, w is the incoherent laser 
pumping, and y is the atom decay rate (we will ignore T2 
dephasing in this example for simplicity). The laser system 
is difficult to solve without approximation since it involves 
both many atoms and large numbers of photons when above 
threshold. Therefore, it constitutes an interesting test-case to 
illustrate the capability of the SU{A) approach. 
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With this in mind, one can now derive a systematic al- 
gorithm for obtaining density matrix elements ^, given 
SU{A) expansion coefficients C™^", ^^,. The procedure is out- 
lined as follows. For each layer of the pyramid [cf. Fig.[Ttc)], 
one may start with a corner element (M and M' maximal) 
and fill out the matrix by successive application of the angu- 
lar momentum lowering operator /_ - c"J (noting that 
pj- = {^++'V+)p) to recursively fill out each row, and 7_p (or 
hermiticity of p) to fill out each column. The layers are filled 
upwards from the base, starting with D"^'"^ n/2 n/2 ~ n/2 
as the corner element of the lowest layer, and finding the cor- 
ner element of higher layers by Gaussian elimination from the 
trace constraint Eq. (|7|. In the supplementary materials, we 
demonstrate explicit application to 3 atoms, with extrapola- 
tion to higher A^ straightforward. 



FIG. 2. (color online). Calculations of laser behaviors described by 
Eq. lEJ with n = 1, 7 = 5, = 1, i/ilToJ = and W = 30. (a) The 
average intracavity photon number (red dots) and spin-spin correla- 
tion (blue squares) as a function of the repumping rate w. The blue 
line is the prediction of the average photon number from the laser the- 
ory [3]. (b) Photon statistics of the laser below threshold w = 4 (blue 
squares) and above threshold w = 8 (red dots), (c) Normalized spec- 
tra of the laser below threshold w = 4 (squares) and above threshold 
w = 8 (dots). The red and green lines are fitted Lorentzian line- 
shapes, (d) Threshold behavior illustrated by the intensity correlation 
g'"*(0) = {a^ aa) I {a^ of- (red dots) and the entropy (blue squares) 
of the whole system as a function of the repumping rate. 

Fig. 12 a) shows the average intracavity photon number of 
the laser as a function of the repumping rate, where the thresh- 
old is evident. This result confirms the conventional laser 
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theory prediction fl- Interestingly, the spin-spin correlation 
{cr^jCr^) above the threshold is directly proportional to the pho- 
ton number, which shows that the collective photon emission 
plays an essential role for the laser In Fig. |2lb), we show that 
the photon statistics of the laser changes from thermal below 
threshold to Poisson above threshold. In Fig.|2|c), we demon- 
strate that the laser linewidth narrows considerably as one 
goes above threshold. Finally, in Fig.|2d), the laser threshold 
behavior is characterized by the intensity correlation g*^'(0) 
and the entropy of the whole system. It can be seen that^<2)(0) 
jumps from two below threshold to one above threshold with 
the entropy increasing and saturating. It is remarkable to have 
an exact solution to this fundamental system that rigorously 
confirms the standard laser theory results. Those results are 
based on various kinds of analytic approximations necessary 
to make the problem tractable. 
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FIG. 3. (color online). Comparison of the second order intensity 
correlation g*^'(0) = {a" a" aa) I (a^ of as a function of the repump 
rate in the steady state superradiance. The green symbols show the 
Monte Carlo results including the statistical errors for A' = 10 atoms 
from Ref. [11]. The red dots show the present calculation using the 
5t/(4) theory And the blue solid line shows the semiclassical results 
from Ref. Illh . Inset: Atom statistics for w = Cir^. The length of 
the bars represents populations of the |5, M) states. 



Ref. Ill III is within error bars. The quantum Monte Carlo sim- 
ulations were of course significantly more numerically inten- 
sive. In the weak pumping limit, the light exhibits strongly 
super-Poissonian fluctuations and deviates remarkably from 
the semiclassical prediction (blue line in Fig. [3ja)). The fail- 
ure of the semiclassical prediction in the weak pumping limit 
indicates that the atoms are in a highly-coiTelated state. To 
reveal the atomic states in this case, we apply the techniques 
of projecting the density operator in the Pq,q^,a-, representation 
onto the 15, M) representation and obtain the populations of 
the atoms on \S,M} states. The inset of Fig.[3]shows explicitly 
that the atoms are mainly pumped into long-lived collective 
subradiant states iH |5 = 0,M = 0) and |5 = 1,M = -1). 
From \S = 0, M = 0), the atoms can only be repumped to 
\S = l,M - 1), from which they rapidly emit two photons 
and relax to IS' - l,M - -1). Therefore, our methods enable 
us to extract detailed information about the atomic states and 
has been essential to reveal the underlying quantum dynamics. 

In conclusion, we have formulated and applied a SU(4) the- 
ory to numerically solve the quantum master equation, which 
has reduced the exponential scaling of the problem to cubic 
in A^. We have developed powerful methods to transform the 
density operator in the SU{4-) basis representation to the \S, M) 
representation. This has enabled us to efficiently diagonal- 
ize the whole density matrix and thus provided complete in- 
formation about the system, including state information and 
functional properties of the density operator We have in- 
cluded lasing and steady-state superradiance as examples in 
order to illustrate the potential for this method. The method 
described here will find numerous applications for simulating 
open quantum systems with large system size. 

We acknowledge stimulating discussions with J. Cooper 
and D. Meiser This work has been supported by the DARPA 
QuASAR program and the NSF. 



As a second example, we apply our approach to steady-state 
superradiance as previously proposed H and demonstrated 
in a recent experiment llioll . The steady-state superradiance 
represents a novel regime of cavity quantum electrodynamics, 
where the highly coherent collective atomic dipole induces an 
extremely nan^ow linewidth for the generated light. The bad- 
cavity mode only plays a role as the source of collective cou- 
pling for the atoms and the definition of the spatial mode for 
the output light [| 1 Ij) . The behavior of this system is also de- 
scribed by a master equation Eq. but in a completely dif- 
ferent parameter regime to the conventional laser For steady- 
state superradiance, the vacuum Rabi splitting is much less 
than the cavity hnewidth, Q. <s: k and equivalently the 
photon number per atom in the cavity is much less than unity. 

We present here calculations of the second order inten- 
sity correlation g^^'(O) in steady-state as a function of the 
repump rate. As shown in Fig. [3ja), the agreement of the 
present calculation and the quantum Monte Carlo result from 
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Supplementary material for 
Simulating Open Quantum Systems using the Simple Lie Group SU(4) 



SU(4) ALGEBRA 



The linear combinations are given by 



In this section, we give explicit expressions for the SU{4) 
superoperators for two-level atoms in terms of the Pauli 
spin-operators. Following Hartmann IJl, 18 superoperators 
0+, O and are defined, where O e {Q, E, M, TV, 1/, 



'U^P 



N ^ N 

N ^ N 

^^P ■= Tj 4 ^ [cr'jp - pcr^) 



N 



^ 1 -CT^ 



N 



N 



N 



(13) 



l-o- 



2^ 2 



2 2 



Although this list, Eq. ( fT3T l. contains 18 operator definitions, 
only 15 of them are independent (it is possible to write A/3, 
1/3, 'V3 in terms of the others). One can demonstrate that 
the 15 remaining superoperators are linear combinations of 
the familiar Gell-Mann matrices that are the generators of the 
SU{4) group, /li , /ii5 In order to see this, consider first 
the fundamental one atom case, where we interpret the 2x2 
density matrix as a 4x 1 vector in the representing vector space 
(i.e. Liouville space). 



a c 
d b 



(14) 



<3+ 



-{Aq ± !/lio), <33 + 



4V3 



^± ^ :^i^6 ± i^i), 13 



1. ^ 



Mi 



± iAs), Mi 



1 



1 



1 



V3 



1 



N± -(/111 ± iMi), Ni -> --/I3 + — -Ai + J -Ais 



4V3 



1/± 



1 



1 



-(.ll + iA2), ^ -.13 



1 , 1 

2 V3 V 6 



(15) 



The commutation relations of the superoperators are given in 
both Ref. ill] and H. We can also identify six SU(2) subalge- 
bras. 



[0^,0 ] ^203, [03,<3±] = ±<3i 



(16) 



so that it is useful to define six corresponding quadratic su- 
peroperators = 0+ + O^ + 0^, which commute with 
Oi. The SU(4) group has 3 Casimir operators, one of which is 
quadratic in the generators, and the others are of higher order 
The quadratic Casimir operator C\ can be expressed in terms 
of superoperators 

Ci = y(0-0++03)+^^j + ^m+2^3f+J(3Q3-21ii-Y3f. 

(17) 

FULLY SYMMETRICAL BASIS FORSt/(4) GROUP 

The fundamental representation of the 5'L'(4) group, 
adapted to serve as basis of the single-spin density matrix, 
is given by u = |1X1|, d = |0)<0|, s = |1><0|, c = |0K1|. 
Higher order representations can then be obtained from the 
fundamental representation and the symmetry type, which is 
described by the Young Tableau. 

The basis for the fully symmetrical case is defined as 



(18) 



which are eigenstates of both O^ and O3 Jill, with eigenvalues 
C?2p(^' = 0(0 + 1)P<'> , O^P'^'^ = 03P<'' , (19) 
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where o 6 {q, cr, m, n, u, v). The eigenvalues are not indepen- 
dent, but can be expressed in terms of a, y,l3, 6: 



q = (a+/3)/2, q^^ {a -(3)12, 

(T = {j + 6)I2, CT^^{y-5)l2, 

m - (a + 5)12, nii - {a - 5)12; 

n = {y+ 13)12, ^ {y - p) 12; 

u = (a + y)/2, Ui - (a - y)/2; 

{6 +13)12, v^^ {6 -13)12. 



(20) 



Then it is straightforward to determine actions of all the rais- 
ing and lowering superoperators on Pq.q^jn, 



Q±Pq,qi,a-j. 

y P 

^l^±P q,q^,(T^ 
'y±Pq,q,,a-j 



{q + qi)Pq,q,±Ua„ 

(cr+cri)Pg 

(m + OT3)^'?±l/2,f/3±l/2,o-3±l/2, 
(n + «3)^'?Tl/2,93±l/2.a-3±l/2, 
(m + U^)Pq±ll2,qi±\il,a-i^ll2, 
(V'+ V3)f^q:l/2,^3±l/2,o-3 + l/2- 



(21) 



We note that the fully symmetrical basis are also eigenstates 
of the quadratic Casimir operator C \ with common eigenvalue 
3MA^-h4)/8. 

Analogous actions for the photon part are the simple har- 
monic oscillator relations: 

a |n) = V« |« ^ 1) ' 

a"'' \n) = Vn + 1 |n + 1 ) . (22) 

|5, Ms){S, M'g I REPRESENTATION 

In order to project the density operator from the SU{4) basis 
onto the \S,M){S,M'\ representation , let us first show that 



M and M' are related to the P, 



'q!q,,iT, by q3 + 0-3 

N _(3) 



M and 



?3 - 0-3 = M'. To see this, defining = 2j=i "j /2, 
could get 



9,93 .o"3 



-ia + y-/3-6)P^\^^,^^ 



(q, + cr^P'i^ 



9.93 .0"3' 



9,93,<r3-^3 ^-ia + 5-p- r)f *1,..3 = (^3 - C^3)f <1,^3' 



and by definition, we have 

J3\S,M}{S,M'\ = M\S,M){S,M'\, 
\S,M}{S,M'\J3 = M'\S,M}{S,M'\. 



(23) 



(24) 



Therefore, the complex coefficients from the P^q^q^^a-^ basis 
contributing to the matrix element for \S ,M){S ,M'\ all satisfy 
qT, + (T-i - M and ^3-0-3 — M' . 

With this in mind, we now describe a systematic algorithm 
to obtain the density matrix elements D"^'"^f^, from the SU{4) 
expansion coefficients C™;^',.^,-,. We illustrate our method by 
considering in detail the elementary case of three atoms. The 



density matrix in the \S , M}{S , M'\ representation is block di- 
agonal in S ; the block matrices for all S can be arranged in 
the shape of a pyramid as shown in Fig. 1(c). For instance, 
the base layer corresponds to 5 = N/2, with the matrix di- 
mension being (A^ + 1)-. The second layer has S = N/2 - I 
and dimension (A^ - 1)-, and so on. Furthermore there are ns 
copies associated with each layer, so that 2s (25 + l)«s - 2^. 
Taking N = 3 for example, there are two layers, S - 3/2 and 
5=1/2 with ^3/2 = 1 and «i/2 = 2, so that the Hilbert space 
dimension is (3 -H 1) H- 2 x (1 H- 1) = 2^ 

The density matrix needs to be built from the bottom layer 
upwards. In the bottom layer, we find that the only element 
contributing to \N/2,N/2}{N/2,N/2\ is f'*v/2A'/2 0- 
left corner is £>;^)2_^/2,iv/2 = ^m,N/2.o- "^^^ ^PP^y low- 
ering operator /_ - o'J to iteratively generate D^'^" ^, 
with M = N/2 - 1, . . . , -N/2. To do this, we need the recur- 
sion relation 



d'^:Ij^,_,^{s,m\p'"-"\s,m' -1)^ 



{S,M\p'"-"L\S,M'} 



V(5 + M'){S - M' + 1) 
{S,M\(^(++^+)p'"'"\S,M') 



y/(S + M')(S - M' + 1) 



(25) 



Therefore, with the actions of the raising and lowering opera- 
tors [Eq. in}], we can derive all D"^'"2n/2M" ^^^^ ^^^^ 
of the bottom layer Using the fact that the density matrix is 



Hermitian and C™^" 



(CTn-^ )*. we could get all the ele- 



9,93,-0-3' 

ments for the first column by Z)^),.^, ;v/2 " (^iv/iro.M')*- 
repeatedly applying the recursion relation [Eq. d25] l1 in each 
row, we then construct the full base layer As an explicit ex- 
ample, we have constructed the bottom layer, i.e. S -3/2 for 
the three atom case. 



13/2,3/2) 
13/2,1/2) 
13/2,-1/2) 
13/2,-3/2) 



(3/2,3/21 <3/2, 1/2| 



'-'3/2,3/2,0 



(3/2,-1/21 <3/2,-3/2| 



V3 



V3 



C 



V3 



V3 



r 



0,0,-3/2 



V3 



c 

(26) 

In order to illustrate the use of the recursion relation, we 
now show how to get D™'" i/t from D'^'^j 3/2- Because 

we have 



V3 



0,0,3/2 



1.-1.1/ 

V3 



3/2,-3/2,0 j 



rw p(s) _ p(s) onH'7/ P*'*' - P*'*' 

•^+^3/2,1/2,0 ~ -^1,1,-1/2 " + ^1/2,1/2,0 ~ -^1,1,-1/2 



^3/2,1/2,1/2 



•C',";.",„„)/V3/V3. 



1,1,-1/2 
'■^3/2,1/2,0 ' ^l/2,l/2,f 

To construct the next layer, we thus find out the top left ma- 
trix element first, and then apply the same procedure as before 
to determine the rest of the matrix elements. Let us first ex- 
amine the three atom case. The 5 = 1/2 layer has two copies, 
each of which is a 2 x 2 matrix. To find the top left element 
D™ 2 J noticing the constraint imposed by the trace of the 
density matrix, we derive 2D"'^2 1/2 1/2 + ^3/2 1/21/2 ~ ^3/2 1/2 
so that Z)'r;^,,/2,i/2 = (2C3';2,i/2,o - Cr/2" i/2.o)/6- By applying 
the same method as in the bottom layer, we construct the block 



7 



matrix for 5 = 1/2 layer 



11/2,1/2) 
11/2,-1/2) 



(1/2,1/21 

' "^3/2.1/2,0 ^1/2.1/2.0 
6 

^].0.-]/2 ^0.0.-1/2 



(1/2,-1/21 



2r"''" -c'"-" 

^3/2.-1/2.0 ^1/2.-1/2.0 



(27) 



Having the top left matrix element for each layer S , we can 
easily construct the (25 + 1) x {2S + 1) block matrix by ap- 
plying the recursion relation based on the angular momentum 
lowering operator Repeated iteration of these steps systemat- 
ically fills in all sites of the pyramid. 



Therefore in general, if we suppose that we have con- 
structed the block matrix for S' > S , the formula to find the 
top left matrix element D™'^' ^ for layer S is 



(28) 



S<S'<N/2 
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